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Abstract 



The identification of primary photons or specifying stringent hmits on the photon 



m 

O 

O ' flux is of major importance for understanding the origin of ultra-high energy (UHE) 

■ cosmic rays. We present a new Monte Carlo program allowing detailed studies of 

Cd ■ 

conversion and cascading of UHE photons in the geomagnetic field. The program 
^ ' named PRESHOWER can be used both as an independent tool or together with 

a shower simulation code. With the stand-alone version of the code it is possible 
to investigate various properties of the particle cascade induced by UHE photons 
interacting in the Earth's magnetic field before entering the Earth's atmosphere. 
Combining this program with an extensive air shower simulation code such as COR- 
SIKA offers the possibility of investigating signatures of photon-initiated showers. 
In particular, features can be studied that help to discern such showers from the 
ones induced by hadrons. As an illustration, calculations for the conditions of the 
southern part of the Pierre Auger Observatory are presented. 
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1 Program Summary 



CPC Program Library Index: 1.1 Cosmic 
Title of program: 

Catalog number: 

Program obtainable: 

Computer on which the program has 
been thoroughly tested: 

Operating system: 

Programming languages used: 

Memory required to execute: 

No. of bits in a word: 

Has the code been vectorized? 

Number of lines in distributed pro- 
gram: 

Other procedures used in 
PRESHOWER: 

Nature of physical problem: 



Rays 

PRESHOWER 1.0 

from P. Homola, e-mail: Pi- 
otr.Homola@ifj.edu.pl 

Intel-Pentium based PC 

Linux, DEC-Unix 

C, FORTRAN 77 

< 100 kB 

32 

no 

1900 

IGRF [1], bessik, ran2 [2] 

simulation of a cascade of particles ini- 
tiated by UHE photon passing through 
the geomagnetic field above the Earth's 
atmosphere. 
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Method of solution: The primary photon is tracked until 

its conversion into e'^e~ pair or un- 
til it reaches the upper atmosphere. 
If conversion occured each individual 
particle in the resultant preshower is 
checked for either bremsstrahlung radi- 
ation (electrons) or secondary gamma 
conversion (photons). The procedure 
ends at the top of atmosphere and the 
shower particle data are saved. 

Restrictions on the complexity of the gamma conversion into particles other 
problem: than electron pair hasn't been taken 

into account 



Typical running time: 



100 preshower events with primary en- 
ergy 10^° eV require a 800 MHz CPU 
time of about 50 min., with 10^^ eV the 
simulation time for 100 events grows up 
to 500 min. 
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2 Introduction 



The existence of ultra-high energy cosmic rays (UHECR) beyond the Greisen- 
Zatsepin-Kuzmin (GZK) cutoff energy [3] is one of the main puzzles in to- 
day's astrophysics. The composition of UHECR, if known, could provide a 
stringent test of models of particle acceleration and propagation towards the 
Earth. Most theoretical efforts are made towards explaining the acceleration 
of protons to ultra-high energies. This is because protons seem to be the 
most promising candidates for the UHECR particles. There are, however, also 
so-called exotic models which postulate decays of supermassive "X-particles". 
A common prediction of these models is that there should be a significant 
photon flux in the ultra-high energy region. Accurate identification of photon 
primaries, measurement of the ultra-high energy photon flux, or specifying the 
upper limit for it, will provide an excellent test for the models of cosmic-ray 
origin. 

Ultra-high energy photons in the presence of the geomagnetic held can convert 
into an electron-positron pair before they enter the atmosphere. Synchrotron 
radiation of this pair gives rise to a cascade of photons. Some of them can 
have energies sufficient for the next conversion. The cascade at the top of the 
atmosphere thus consists mainly of photons and two or more electrons. We 
will call this cascade a "preshower" since it originates and develops above the 
upper atmosphere, i.e. before the "ordinary" shower development in the air. 

Particles of a preshower induce subshowers after entering the atmosphere. A 
superposition of these subshowers will be seen by detectors as one extensive 
air shower, whose profile is expected to differ from the profile of a shower 
induced by a single unconverted photon. 
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Forthcoming experiments with large collecting apertures (Pierre Auger Ob- 
servatory [4], EUSO [5]), that are designed to detect cosmic rays of ultra-high 
energies, give a good opportunity to test these models. A photon-induced 
shower might have observable features that allow to distinguish it from show- 
ers initiated by hadrons. 

In this paper we describe the PRESHOWER program that allows extensive 
studies of properties of photon-induced showers. The program performs a sim- 
ulation of the propagation of photons before they enter the Earth's atmosphere 
and can be used both as an independent analysis tool as well as a part of an 
air shower simulation package (e.g. CORSIKA [6]). 

With the PRESHOWER program working in stand-alone mode we performed, 
as an illustration, an analysis of basic properties of preshowers for the magnetic 
conditions of Malargiie in Argentina - the site of the Southern Pierre Auger 
Observatory [4]. The results of these studies, which are presented in Sec. 5, 
are in general agreement with other analyses [7,8,9,10,11,12,13]. 

By linking our program to an air shower simulation code one is able to 
study showers initiated by various primaries. We checked this by combining 
our PRESHOWER code with the CORSIKA air shower simulation package 
which simulates air showers initiated by the preshower particles. We simu- 
lated photon-induced showers of different energies and different directions at 
the Southern Pierre Auger Observatory. Although the interaction in the ge- 
omagnetic field changes significantly the shower properties so that primary 
photon events look more similar to hadronic showers, a discrimination be- 
tween primary photons and hadrons still seems possible. 

The plan of the paper is the following: The relevant effects and the formulas 
applied in the simulation are summarized in Section 3. The structure of the 
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program is discussed in Section 4. As an illustration, in Section 5 results 
obtained with the code for the case of the Auger Observatory are discussed. 
The paper is concluded in Section 6. 

3 Preshower physics 

At energies above 10 EeV, in the presence of the geomagnetic field, a photon 
can convert into an e^e~ pair before entering the atmosphere. The resultant 
electrons will subsequently lose their energy by magnetic bremsstrahlung. The 
emitted photons can convert again if their energy is high enough. In this way, 
instead of one high energy photon at the top of atmosphere, a number of less 
energetic particles, mainly photons and a few electrons, will emerge consti- 
tuting a preshower. A superposition of subshowers induced by these particles 
should be seen by detectors as one extensive air shower which reaches its max- 
imum earlier than a shower induced by a single photon of the same primary 
energy. This is also due to the Landau-Pomeranchuk-Migdal (LPM) effect [14] 
which is stronger for photons of higher energies. Both photon conversion and 
magnetic bremsstrahlung are strongly dependent on the transverse component 
of the geomagnetic field. In the following, we concentrate on these effects, give 
the main formulas and discuss their application. 

3.1 Geomagnetic field model 

Since the geomagnetic field is responsible both for gamma conversion and 
for bremsstrahlung of electrons, it is very important to model it with ade- 
quate accuracy. We use the International Geomagnetic Reference Field model 
(IGRF) [15] for n — lQ,n being the highest order of spherical harmonics in the 
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scalar potential. To calculate the field according to the IGRF model we use the 
numerical procedures written by Tsyganenko [1]. When choosing n — 1, the 
output of these procedures is equivalent to a dipole model with |M| = 8 • 10^^ 
G-cm^ and with poles at 79.3°N, 71.5°W and 79.3°S, 108. 5°E [15]. 



3.2 Magnetic pair production: 7 — e+e 



The number of pairs created by a high-energy photon in the presence of a 
magnetic field per path length dr can be expressed in terms of the attenuation 
coefficient a{x) [16]: 

ripairs = riphotons{i " exp[-Q;(x)cir] } (1) 



where 

a{x) = OMaemmec/h){B^/B,r)T{x) (2) 



with being the fine structure constant, x = 0.5(/ii//meC^)(Sx/Scr), B± 
is the magnetic field component transverse to the direction of the photon's 
motion, Bcr = mic'/eh = 4.414 x 10^"^ G and T(x) is the magnetic pair 
production function. T(x) can be well approximated by: 



T{x)^0.16x-'K\/,{—), 



(3) 



where Ki/^ is the modified Bessel function of order 1/3. For small or large 
arguments T{x) can be approximated by 



T{X) = 



0.46exp(-|^), x«l; 



(4) 



0.60X 



-1/3 



X> 1- 
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Representative values of T(x) are listed in Table 1. 

We use Eq. (1) to calculate the probability of 7 conversion over a small path 
length dr: 

Pconv{r) = 1 - ex.p[-a{x{r))dr] ~ a{x{r))dr (5) 
and the probability of 7 conversion over a longer distance R: 

R 

Pconv{R) = 1 - exp[- j a{x{r))dr] (6) 



3.3 Magnetic bremsstrahlung 

After photon conversion, the electron-positron pair propagates. The energy 
distribution in an e'^e~ pair is computed according to Ref. [17]: 

da{£, x) aemiriecB^ 3^/2 [2 + e{l - e)] 



de hBcr 97rx 6{1 — e) 



-K2/3 



3xe{l - e) 



(7) 



where s denotes the fractional energy of an electron and the other symbols 
were explained in the previous chapter. The probability of asymmetric energy 
partition grows with the primary photon energy and with the magnetic field. 
Beginning from x > 10, the asymmetric energy partition is even more favored 
than the symmetric one. 

Electrons traveling at relativistic speeds in the presence of a magnetic field 
emit bremsstrahlung radiation (synchrotron radiation). For electron energies 
E ^ rrieC^ and for B± <^ Bcr, the spectral distribution of radiated energy is 
given in Ref. [18]: 



9 



where ^ = {3/2){B±/ Bcr){E/'meC^), E and rUe are electron initial energy and 
rest mass respectively, ^5/3 and ^2/3 are modified Bessel functions, and y is 
related to the emitted photon energy hu by 

/;// dy E 

' ^ C{E - hv) ' dijw) ^ i{E - hvf 

The total energy emitted per unit distance is (in CGS units) 

2 °° 



with ro being the classical electron radius. For our purposes we use the spectral 
distribution of radiated energy defined as 

where dN is the number of photons with energy between hv and hv + dihv) 
emitted over a distance dx. From Eqs. (8), (9), (10), and (11) we get ^ 

nB,,EM) = lrlBl{-^p(h.)^^^^ . (12) 



Provided dx is small enough, dN can be interpreted as a probability of emitting 
a photon of energy between hv and hv + d{hv) by an electron of energy E 
over a distance dx. In our simulations we use a small step size of = 1 km. 
The total probability of emitting a photon in step dx can then be written as 

Pbrem{B±, E, hv, dx) = j dN = dx j I{B^, E, H^^^ ■ (13) 





^ Expression (12), valid for all values of hu, is equivalent to Eq. (2.5a) in Ref. [16]. 
A simplified form of distribution (12) is given by Eq. (2.10) in Ref. [16], however it 
can be used only for hv <C E. 
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The shape of the bremsstrahlung spectrum for two different energies is shown 
in Fig. 1. The three curves in each plot correspond to different values of the 
transverse magnetic field: 0.3 G, 0.1 G and 0.03 G from the uppermost to the 
lowest curve respectively. The energy of the emitted photon is simulated ac- 
cording to the probability density distribution dN/d{hi') obtained from Eq. 13. 



3.4 Related physics topics 



In this section we shortly discuss the importance of other physics effects con- 
nected to our main subject. These are: 

• The geomagnetic field defiects the electrons from the initial preshower di- 
rection. To evaluate this defiection we use the equation 

— — = evB± ; rrirei = (14) 

where v is the speed of electron and e is the elementary charge, to find the 
radius of curvature R. If it! is much larger than the electron path length L 
and assuming constant values of S^, E and v ~ c, the linear displacement 
from the preshower core in the plane perpendicular to the initial direction 
is well approximated by: 

Ax ^ — = ^ — 15 

2R E ^ ^ 

For typical values like B± — 0.1 G, L — 1000 km (see Section 5) and 
relatively small energy of E — 10^^ eV, we get R ~ 10^^ km which gives 
Ax = 0.5 cm. Assuming a larger value for L, one should keep in mind that 
at higher altitudes B is lower and E is higher. Allowing for other possible 
combinations of B±, E and L one comes to a very safe conclusion that in 
any case Ax < 1 m. 
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The initial angular spread of particles in both pair production and brems- 
strahlung can be approximated by ~ vrtf.jE [19]. Even for energies as 
low as 10^^ eV one has A^^ 10~^^ which even for extremely long paths like 
for instance 10000 km results in negligible linear deviations at the top of 
the atmosphere of Ax ~ 10~^cm. 

A 7-ray traveling through the vicinity of the Sun can cascade in the mag- 
netosphere of the Sun where the field is much stronger than near the Earth. 
The particles originated in conversion near the Sun can be significantly de- 
fiected from their initial direction by the solar magnetic field. The cascade 
of such particles after reaching the Earth's atmosphere will induce an air 
shower of large lateral extent which might be well distinguishable from pro- 
ton or iron showers. However, the predicted rate of such events for primary 
energies of order 10^^ eV is very low: about 1 event per 10 years for the Pierre 
Auger Observatory [20]. In our simulations we do not consider cascading in 
the Sun's magnetosphere. 

The magnetic field embedded in the solar wind is of order 10"^ G [21] and it 
influences the Earth's magnetosphere at altitudes of order 8-10/?e. Even for 
a primary energy as high as 10^^ eV gamma conversion takes place within 
2i?0 above the sea level (see Section 5) where the Earth magnetic field is 
not less than 10~^ G. Thus, the solar wind can be neglected in our studies. 
The preshower electrons traveling at speeds v < c are delayed with respect 
to the photons. This delay however is very small and can be neglected. For 
an extreme case of a 10^^ eV electron traveling a path of 20000 km, the 
delay is of order 10~^^ s. The delay of the particles not parallel to the main 
direction is also negligible (see the above discussion on the lateral spread of 
a preshower). 
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4 Structure of the program 

The following subsections refer to the stand-alone mode of the PRESHOWER 
program. The PRESHOWER-CORSIKA mode requires only the main proce- 
dure generating events, preshw, which is described in subsection 4.3. All the 
input data in the PRESHOWER-CORSIKA mode are provided via the main 
CORSIKA input file, the standard output is printed within the CORSIKA 
output, and the multirun.dat file is not returned, since there is only one run 
per one primary trajectory allowed in CORSIKA. There is also a different ran- 
dom number generator used in the PRESHOWER-CORSIKA mode than in 
the stand-alone one. For more information on the PRESHOWER-CORSIKA 
mode the user can refer to the latest CORSIKA manual [22] or contact the 
corresponding author. 

4-1 The source code and other files in the PRESHOWER packet 
The PRESHOWER code consists of the following files: 

preshw. c contains the main procedure generating preshowers, preshw, and 
the other procedures in C that are called by preshw. 

prog.c is used only in the stand-alone mode. It reads the input file INPUT 
and calls the preshower procedure preshw. 

igrf.f is an external procedure for calculation of the geomagnetic field ac- 
cording to the IGRF model (author: N.A.Tsyganenko [1]), it is called by 
preshw. 

rmmard.f is a dummy Fortran routine used only in the stand-alone version 
for the proper linking of the program. In the PRESHOWER-CORSIKA 
mode the actual rmmard is responsible for random number generation. 
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For the user's convenience there are some more files included in the PRESHOWER 
packet: 

Makefile compilation and linking commands under linux. 

README.txt short information on the content of the PRESHOWER packet 

and on the installation process. 
INPUT typical input data read by the program working in the stand-alone 

mode. In the PRESHOWER-CORSIKA mode all the input data required 

by preshw are read directly from the main CORSIKA input file. 
part_out.dat the example of a detailed output file containing information 

on the preshower particle data being the result of running the program with 

the input INPUT. 

multirun.dat the example of an output file summarizing the program run 
with the input INPUT; this file is produced only in the stand-alone mode. 

out.txt the standard output produced by the program in the stand-alone 
mode with the use of the input data from INPUT. 

4-2 Input and output data 

The input data are read by the main program from the file INPUT and then 
passed to the preshw procedure. These data are 

(1) the primary photon energy, [GeV] 

(2) the zenith angle of the shower trajectory in the local reference frame (see 
Sec. 5), [deg] 

(3) the azimuth angle of the shower trajectory in the local frame, [deg] 

(4) the top of the atmosphere - e.g. the value obtained from CORSIKA, [km 
a.s.l.] 

(5) geographical position of the observatory: longitude [deg] 
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(6) geographical position of the observatory: latitude [deg] 

(7) the year of observation (needed for calculation of the field) 

(8) number of runs 

An exemplary set of the input data can be found in the file INPUT. 

There are two output files produced by the PRESHOWER program: part_out . dat 
and multirun.dat. Both files begin with a header containing all the crucial 
run information, then comes the description of columns and then the data. The 
data in part_out.dat are presented in two columns: particle id (1 - photons, 
2 - positrons, 3 - electrons) and particle energy at the top of the atmosphere 
(in eV). The columns in multirun.dat are explained as follows: 

(1) the run number 

(2) altitude of the first photon conversion in km a.s.l. 

(3) total number of particles at the top of the atmosphere 

(4) number of photons 

(5) number of electrons 

(6) maximum electron energy 

(7) maximum photon energy 

(8) fraction of energy carried by electrons 

(9) total energy carried by preshower particles = primary photon energy 

(10) total energy carried by electrons 

(11) total energy carried by photons 

The information printed as standard output begins with the same header as 
the output files and is followed by the specific information on each run (all 
energies given in eV) : the run number, the altitude of each conversion together 
with the energy partition information, preshower summary at the top of the 
atmosphere: total energy carried by photons (Etot_g), total energy carried by 
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electrons (Etot_e), maximum photon energy (Eg_max), maximum electron 
energy (Ee_max), number of particles (n_part), number of photons (n_phot), 
number of positrons (n_e+) and number of electrons (n_e-). The example of 
a standard output is the file out.txt included in the PRESHOWER packet. 

4-3 Generation of events 

After reading the input parameters (see Sec. 4. 2) the event generation proce- 
dure preshw is called, within which all the calculations are performed. In the 
initial stage all the necessary parameter transformations are applied and the 
preshower trajectory vector is computed. The propagation path length is set 
to be 5i?® long, which e.g. for vertical showers means the simulation starting 
altitude of 5Rq above sea level. All the data concerning the particles in the 
preshower are stored in the array part_out which at the beginning has only 
one entry: the primary photon. In steps of 10 km the transverse magnetic field 
is computed by the procedure B_transverse, then the probability of conver- 
sion is found using Eq. (5). The step size has been checked to be sufficiently 
small. If the conversion takes place, the particle array is updated by replacing 
the photon with an e"*" and adding an e~ at the array end. The energy partition 
fraction is found by the function ppfrac. To assure accurate bremsstrahlung 
simulation, beginning from the first gamma conversion the step of 1 km is 
used for all particle types. In every step, for each photon in the array, we 
check again for the conversion effect. In case of electrons the bremsstrahlung 
radiation is simulated using function brem. The probability of emitting a pho- 
ton is calculated from Eq. (13) with the integration starting from hu = 10^^ 
eV and logarithmic steps of 0.01. Bremsstrahlung photons of energies lower 
than 10^^ eV have a very small infiuence on the air shower evolution (see also 
Section 5) and hence they are neglected. 
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Once a bremsstrahlung photon is emitted, its energy is determined by the 
probability distribution dN/d{hv) obtained again from Eq. (13), and the en- 
ergy of the radiating electron is diminished correspondingly. Simulations are 
finished when the (adjustable) altitude of the top of the atmosphere is reached 
and all preshower particle parameters are written out (in the stand-alone 
mode) or passed to CORSIKA. 



4-4 Compilation, linking and running the program 

In order to compile, link and run PRESHOWER under Linux one should first 
create a file nr_f un . c in the PRESHOWER home directory and place the 
header: 

#include <stdio.h> 
#include <math.h> 

at the beginning of this file. Then, below the header, the source code of the 
Numerical Recipes [2] procedures ran2 and bessik (together with the proce- 
dures called by bessik: beschb, chebev and nrerror) should be pasted. To 
compile and link the program one just types "make" in the PRESHOWER 
home directory. For Unix compilation and linking one needs to replace the 
proper commands in Makefile. Once the data in the INPUT file are appropri- 
ate, the program can be started by typing "./preshower". The compilation and 
linking process can be verified by running the program with the exemplary 
INPUT file included in the packet. If the program runs correctly, it should re- 
turn the same output files as those included in the packet (part_out.dat, 
multirun.dat) and the standard output as in out.txt. 
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5 Illustration of preshower simulations 

The results presented in the following are obtained for the magnetic conditions 
of the Southern Auger Observatory in Malargiie, Argentina (35.2°S, 69.2°W). 
Other geographical positions could easily be adopted. The shower trajecto- 
ries are given in the local frame where the azimuth increases in the counter- 
clockwise direction and (p — 0° means a shower coming from the geographical 
North. 

To visualize how the transverse magnetic field {B±) can change along the 
shower trajectory, we have chosen some representative directions along which 
a particle encounters B± of significantly different strengths: strong, medium 
and weak. The "strong field direction" is the one along which B± is maximal 
compared to other directions with zenith angles 9 < 60° and is defined by 
the angles 9 = 60° and (p = 177°. Similarly, the weak B± trajectory ("weak 
field direction") was chosen to resemble the local field direction at Malargiie 
of ^ = 50°, (j) — 357° [15]. The "medium field direction" is a vertical trajectory 
{9 — 0°). In the following we will refer to these directions. 

5.1 Geomagnetic cascading 

In Fig. 2 we show how B± changes with altitude for the strong, medium and 
weak field directions. The trajectories end at the chosen top of the atmosphere. 
In addition to the altitude dependence of B±, large directional differences are 
visible. Of course we expect the strong field direction to coincide with the 
direction of a stronger preshower effect. 

In Fig. 3, the altitudes of first photon conversion are shown for different ener- 
gies in the strong field direction. The percentage values shown for each curve 
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denote the total conversion probability up to the top of atmosphere, which 
is adjusted here to the starting altitude of CORSIKA at 112 km. It can be 
seen that for the magnetic field above Malargiie, the pair production effect 
becomes important for energies above 5 x 10^^ eV. For all energies of interest, 
even in case of photons at 10^^ eV, the conversion altitude is not higher than 
13000 km, i.e. about 2R^ a.s.l. Thus, the range of interest of the magnetic 
field strength for the case of Malargiie is 0.3 G -^ 0.01 G (see Fig. 2). 

The distribution plotted in Fig. 3 is consistent with the results of Ref. [9]. 
The integrated conversion probability computed with Eq. (6) for different 
7 energies and for the two trajectories of strong and weak field directions is 
plotted in Fig. 4. For photons of energies between 5 x 10^^ eV and 5 x 10^° eV a 
change in the arrival direction, i.e. a change in the magnetic field component 
perpendicular to the photon trajectory, may imply a dramatic increase or 
decrease of the conversion probability. Thus, for these energies we expect a 
strong directional dependence of the preshower effect which might serve as 
characteristic fingerprint of photons as primary cosmic rays. 

The directional dependence of the gamma conversion probability is shown in 
more detail in Fig. 5 for different arrival directions and primary energies. The 
conversion effect starts to be important for directions towards the magnetic 
South pole at photon energies around 50 EeV. For energies of about 150 EeV 
the conversion is very probable for almost every direction. The photons com- 
ing from azimuthal angles around 150° at large zenith angles convert more 
likely than the others. To explain this, one notes that for a given altitude the 
magnetic field is strongest around the geomagnetic axis and of course the field 
becomes stronger with decreasing altitude. Since the photons coming from the 
geomagnetic South, which corresponds to the azimuth about 150° in case of 
Malargiie, cross the very vicinity of the geomagnetic axis, they must encounter 
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the stronger magnetic field than the others and hence their conversion prob- 
ability is higher. The larger zenith angle for these photons means that the 
strong field region is crossed at the lower altitudes where the field is stronger 
and the conversion process is more effective. The plots in Fig. 5 are consistent 
with the values given in Ref. [13]. 

All the presliowers discussed here are presented in their final stage of develop- 
ment at the "top of the atmosphere" which in our simulations is assumed to be 
112 km a.s.l. - where in the CORSIKA code shower simulations are started. 

In Fig. 6 mean energy spectra are plotted for 1000 preshowers initiated by 
photons of energy 10^0 eV arriving at the top of atmosphere from the strong 
field direction for Malargiie. In the top row we show the distribution of photons 
(left) and the distribution of photons weighted by their energies (right) . We 
start the simulations of the bremsstrahlung radiation for the energy threshold 
of 10^^ eV so the photons of lower energies are absent in the distribution. In 
this way we take into account the majority of the photons even if those of 
energies lower than 10^*^ eV may have little impact on the longitudinal profile 
of the resultant air shower. The total energy fraction carried by such photons 
is very low (see the top right plot in Fig. 6) and the induced subshowers fade 
very quickly. In the simulation run presented in Fig. 6, only 78 out of 1000 
primary photons did not convert into e"'"e~, these are seen in the top right 
plot as the narrow peak at E — lO^'^ eV. 

The distributions of the electrons in the same sample of 1000 preshowers 
are shown in the bottom row of Fig. 6. Again, the spectrum of particles is 
plotted to the left and the spectrum weighted by energy to the right. Only in 
three cases out of 1000 a bremsstrahlung photon converted into e'^e~ so the 
total number of electrons is 1850 and not 1844. When the conversion happens 
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at a lower altitude, the electrons radiate less bremsstrahlung before entering 
the atmosphere. These cases are seen in the spectrum as a tail towards high 
energies. 

Some characteristics of the preshowers shown in Figure 6 are summarized in 
Table 2 together with some other exemplary preshower cases. The conversion 
probabilities are consistent with those in Figure 5. For higher initial energies 
gamma conversion takes place at higher altitudes (see Fig. 3) and the resultant 
electrons radiate more bremsstrahlung photons so the number of particles is 
larger. Additionally for energies as high as 10^^ eV there is more than one 
gamma conversion per each preshower which also enlarges Npart- Comparing 
strong Bj^ and weak Bj^ directions for a given energy, the large directional 
dependence of preshower evolution can be noted. In case of Eq — 10^° eV from 
the weak Bj^ direction we expect unconverted photons initiating cascades with 
a large depth of shower maximum X^nax due to the LPM effect. The strong 
B±^ forces the initial energy to be split into about 500 particles of which the 
highest energy particle typically is a photon with the energy about 5 times 
lower than the initial one. This means that the longitudinal profile of the 
resultant air shower will have the maximum higher in the atmosphere, which 
will make it more similar to a proton-like profile. The energy splitting in this 
case will compete with the LPM effect. This splitting effect for Eq = 10^^ eV 
for the strong B± direction is even more dramatic: the highest energy particle 
in the preshower is more than 16 times less energetic than the primary photon. 
The maximum 7 energy in the preshower for the weak field direction is higher 
than for the strong field since the conversion probability is lower in a weaker 
field (see Fig. 4). Detailed distributions of <Npart>, <A^e+e->; <max E^> 
and <max Ee> for Eq — 10^^ eV and the strong B± are shown in the top left, 
top right, bottom left and bottom right plots of Figure 7, respectively. 
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As mentioned above, the number of particles in the preshower is connected 
to the altitude of the first gamma conversion. This correlation is presented 
in Figure 8. The higher altitude of conversion gives more particles at the top 
of the atmosphere because the electron-positron pair radiates bremsstrahlung 
over a longer path. The increase in particle number is slightly larger at low 
altitudes because of the stronger magnetic field. 

For the same reasons the total energy fraction carried by the electrons should 
be smaller for higher conversion altitudes. This dependence is shown in Fig- 
ure 9. One notes that for high altitudes of the first conversion only about 1% 
of the initial energy of 10^° eV remains in the electrons at the top of the atmo- 
sphere, while all the rest is radiated into photons. With decreasing altitudes 
the energy radiated into photons decreases. The three points in excess of the 
general trend in Figure 9 are the cases where the primary photon converted 
at high altitude and one of the resulting bremsstrahlung photons converted 
close to the top of the atmosphere - at 228, 442 and 133 km a.s.l respectively 
for the points from right to left. The threshold for 7 conversion in this case is 
about 5 • 10^^ eV, so at least one member of e^e~ pair had an initial energy 
greater than 10^^ eV, and since its birth altitude was relatively low, there was 
little time for radiation. 

5.2 Atmospheric cascading: PRESHOWER and CORSIKA 

The results presented above give a general idea about the characteristics of 
preshowers initiated by UHE photons before entering the atmosphere. The 
next step is to simulate air showers induced by different preshowers or un- 
converted photons and compare the results with the proton-induced showers 
in order to investigate the signatures of UHE photons that will be measur- 
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able by the Pierre Auger Observatory. For this purpose we connected the 
PRESHOWER program to CORSIKA, a commonly used air shower simula- 
tion code. 

The connection between the two codes is organized as follows. First, the dis- 
tribution of the preshower particles at the top of the atmosphere is simulated 
with the original preshower code, then the output is passed to CORSIKA and 
finally the resulting air shower is simulated as a superposition of the subshow- 
ers initiated by the preshower particles. 

Electromagnetic interactions are simulated in CORSIKA using an adapted 
version [23] of the EGS4 code [24], which includes the Landau-Pomeranchuk- 
Migdal effect [14]. To reduce the computational effort in CPU time, the tech- 
nique of particle thinning [25], including weight limitation [26,27], is applied. 
More specifically, a thinning level of 10^^ has been chosen, and particle weights 
are limited to 10~^£^o/GeV, with Eq being the primary energy. This keeps 
artificial fiuctuations introduced by the contribution of individual particles 
sufficiently small for the analysis of shower maximum depths [28]. 

In Table 7 we compare X^ax and RMS of X^ax for gamma induced showers of 
different primary energies and arrival directions. The X^^ax of proton-induced 
showers using the high-energy hadronic interaction model QGSJET 01 [29] 
is 820 (870) g/cm^ for lO^o (lO^^) eV, and using the SIBYLL 2.1 model 
[30] 855 (915) g/cm^. The RMS in all cases is about 50-60 g/cm^. Both 
for Eq — 10^^'^ eV in strong B± direction and for Eq — 10^° eV in weak 
B± direction, almost every photon remains unconverted when entering the 
atmosphere, resulting in deep (X^nax) and large fluctuations due to the LPM 
effect. A comparison to the longitudinal profile of hadronic primaries shows 
that unconverted primary photons might be well distinguishable from primary 
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protons or nuclei even on event-by-event basis. 

Almost all the photons at lO^'^ eV arriving from the strong field direction 
convert into e'^e~ pairs. As the most energetic bremsstrahlung photons have 
energies several times lower than the primary photon energy, the subshowers 
develop faster. On average, the showers reach their maxima earlier than the 
one induced by an unconverted photon, but still not as early as hadronic 
showers. For larger statistics of 10^° eV events and if a substantial fraction 
of UHE cosmic rays are photons, one expects to see, as described above, the 
directional anisotropy in {X^ax) and 'RM.^{X^ax) and on this basis one may 
draw conclusions on the primary composition. The absence of this dependence 
would allow us to set upper limits on the photon flux. 

At the primary energy of 10^^ eV all photons convert, whatever the arrival 
direction. We still see the directional dependence of {X^^ax) but it is not as 
pronounced as previously. The fluctuations of X^ax in this case are signifl- 
cantly lower (by about a factor 4) than for 10^° eV primaries. Table 7 shows 
that at 10^^ eV it is more difficult to distinguish a photon primary on the basis 
of the Xmax value from a proton one on an event-by-event basis. It is interest- 
ing, however, to note the directional dependence of the elongation rate. From 
Table 7 we see that at Malargiie, for the strong direction, the elongation 
rate of photon-induced showers between 10^° eV and 10^^ eV is much less than 
for proton or iron showers. For the weak Bj^ we even have a negative elongation 
rate. This is because the preshowering effect for photons at 10^^ eV splits the 
initial energy into energies less than 10^° eV and at this energy level, for the 
weak Bj^ direction, almost all the primary photons remain unconverted and 
they induce air showers with deeper Xjnax- Low or negative elongation rates 
could be an additional good signature of photon showers. 
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More detailed simulations of the signatures mentioned above will allow us to 
quantify the expected sensitivity of the Auger experiment to a photon fraction 
in UHE cosmic ray flux. Also other observables will be studied, especially the 
muon content and the features of the lateral distribution measured by the 
surface detector, since for the surface measurements the event sample will be 
about 10 times more numerous than for the fluorescence detector. 

6 Summary 

We presented a method for studying the signatures of UHECR gamma show- 
ers and the practical implementation of this method in the PRESHOWER 
program. This program was tested extensively both in stand-alone mode and 
as a part of CORSIKA and was checked for consistency with results quoted 
in other publications. The preshower option will be accessible in the future 
release of CORSIKA. 

With the code working in stand-alone mode, we studied the preshower char- 
acteristics for different primary energies and arrival directions. The strong 
directional and energy dependence of the preshower effect requires a concise 
treatment in the simulations. Anisotropics in the preshower cascade translate 
into anisotropics of air shower observables such as {Xmax) and KW^i^Xj^ax) . 
First results obtained with PRESHOWER combined with the CORSIKA air 
shower simulation code for the conditions of the Auger Observatory confirm 
the expectation that photon showers can be distinguished from hadron ones. 
Further analyses including other observables such as the muon shower content 
are in progress. 

The new experiments with large collecting apertures give promising prospects 
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for acquisition of a substantial amount of extensive air shower data in the 
UHE range of cosmic-ray spectrum. An important step towards explaining 
the origin of UHE cosmic rays will be the observation of primary UHE pho- 
tons or specifying stringent upper limits for the photon flux. This requires a 
detailed understanding of the UHE photon shower characteristics. The pre- 
sented PRESHOWER code can serve as tool for such simulations. 
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A Compilation of procedures/functions 



preshw(*id, *E_gamma, *the_loc, *phi_loc, *toa, *glong, *glat, 
igrf_year, *iiprint, *corsika, *nruns, part_out [10000] [7] , 
r_zero, *r_fpp, *N_part_out) : 

This is the main procedure of PRESHOWER which simulates the prop- 
agation of an UHE photon before its entering the Earth's atmosphere. It 
takes into account the effect of 7 conversion into e^l~ pair and subsequent 
bremsstrahlung radiation of the electrons in the geomagnetic field. As a re- 
sult a bunch (a preshower) of particles is obtained, mainly photons and a 
few electrons, instead of the primary photon. The information about all the 
particles in the preshower is saved (in the stand-alone mode) or returned 
to CORSIKA (PRESHOWER-CORSIKA mode). The output parameters 
r_zero, r_f pp and N_part_out are required by CORSIKA only - they are 
not used by the stand-alone version. 
INPUT: 

(1) *id: primary particle type: always id=l (a photon) 

(2) E_gamma: primary energy [GeV] 

(3) *the_loc, *phi_loc: zenith and azimuth angles (in radians) of the 
shower trajectory given in the local coordinate system 

(4) *toa: the top of the atmosphere [cm] 

(5) *glong: longitude position of experiment [deg] (Greenwich glong = 0°, 
eastward is positive) 

(6) *glat: latitude position of experiment [deg] (North Pole: glat=90°. 
South Pole: glat=-90°) 

(7) *igrf_year: the year in which the magnetic field is to be computed - 
this parameter is an input for the igrf procedure 

(8) *iiprint: print flag (adjustable only in PRESHOWER-CORSIKA, in 
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the stand-alone version iiprint=l) 
= 0: suppress printing from preshw; 
= 1: enable printing from preshw 
(9) *corsika: mode flag 

= 1: PRESHOWER-CORSIKA: random generator from CORSIKA 
is used (rndm); 

= 0: stand-alone version: random generator from Numerical Recipes 
(ran2) is used 
(10) *nruns: number of runs 
OUTPUT: 

(1) part_out [10000] [7] : the array to store the preshower particle data, 
maximum number of particles: 10000 

part_out[n][0]: particle type: 1 - photon, 2 - e"*", 3 e" 
part_out[n][l]: particle energy [GeV] 
part_out[n][2-6]: not used 

(2) *r_zero: height above sea level of the simulation starting point given 
in [km] above sea level 

(3) *r_fpp: height of the first pair production (0 for surviving photons) 
given in [km] above sea level 

(4) *N_part_out: number of preshower particles at the top of the atmo- 
sphere (number of entries in part_out array) 

bessik(x, xnu, *ri, *rk, *rip, *rkp): 
for argument x computes modified Bessel function rk of fractional order 
xnu; [2] 

getrandCseed, mode): 
for a given seed seed returns a random number from interval (0,1) choosing 
an appropriate generator according to the mode flag: ran2 for the stand- 
alone mode (mode=0) and rmmard for the CORSIKA-PRESHOWER (mode=l) 
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float kappa (x): 

an auxiliary function used in the bremsstrahlung procedure brem 
brem(B_tr, E_e, E_sg): 

computes bremsstrahlung radiation for the transverse magnetic field B_tr, 

electron energy E_e and bremsstrahlung photon energy E_sg 
B_transverse (B , tn): 

getting the geomagnetic field B component which is perpendicular to the 

preshower trajectory tn 
ran2(seed) : 

the random number generator used only in the stand-alone mode; [2] 
cross(a, b, c): 

returns cross product c of vectors a and b 
dot (a, b): 

returns dot product of vectors a and b 
norm(a, n): 

the input vector a is normalized and returned as vector n 
normv(a) : 

returns the length of vector a 
glob21oc(glob, loc, theta, phi): 

converts global cartesian coordinates glob to the local ones (loc) for the 

geographical position theta (colatitude) and phi (longitude) 
sph2car(sph, car): 

converts spherical coordinates sph to the cartesian ones car. 
car2sph(car, sph): 

inverse to sph2car 
locsphC2glob(locsph, glob, theta, phi, sing, cosgm, sitec): 

converts local spherical coordinates locsph (CORSIKA frame of reference) 

to the global cartesian coordinates glob, using geographical coordinates 
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of the observatory (theta - colatitude, phi - longitude), global cartesian 
coordinates of the observatory (sitec) and sine and cosine of the declination 
angle (sing, cosg) 

pppCefrac, B_tr, E_gO): 
this auxiliary function is used within ppf rac 

ppfracCs, mode, B_tr, E_gO): 
determines the energy partition after gamma conversion, the computations 
are performed according to Ref. [17]; s - seed, mode - determines the type 
of the random number generator, B_tr - transverse magnetic field, E_gO - 
primary photon energy 

igrf_(year, n, r, theta, phi, Br, Btheta, Bphi): 
external Fortran procedure written by Tsyganenko [1]; computes magnetic 
field according to the IGRF model; used within b_igrf ; year - the year of 
the experiment, n - the highest order of the spherical harmonics in the scalar 
potential, r, theta, phi - spherical coordinates of the point at which the 
field is to be found, Br, Btheta, Bphi - spherical coordinates of the mag- 
netic field at given position 

b_igrf(year, r, theta, phi, bear): 
for a given year and position (r, theta, phi) computes the magnetic field 
bear in global cartesian coordinates according to the IGRF model 

bsph2car (colat , longit, bsph, bear): 
converts the magnetic field from spherical coordinates (bsph) into cartesian 
ones bear; used only in b_igrf 

rndm ( ) : 

calls CORSIKA random number generator, used only in the PRESHOWER- 
CORSIKA mode 
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Table 1 

The magnetic pair production function T{x)- The values below differ from the ones 
given in Ref. [16]. The computation of the limiting cases of (Eq. (4)) or checking 
the K function values in mathematical tables proves that the numbers given here 
are correct. 



Y 


0.2 


0.4 


1.0 


•5.0 


10.0 


.30.0 


100.0 


T{X) 


5.9e-4 


1.6e-2 


9.9e-2 


2.2e-l 


2.1e-l 


1.9e-l 


1.3e-l 



Table 2 

Exemplary preshowers for Malargiie magnetic conditions. For each type of the 
preshower determined by the initial photon energy Eq and the arrival direction, in 
columns 3-7 respectively shown are: the fraction of converting primary photons and 
the average values (± RMS) for the number of particles, the number of electrons, the 
maximum photon energy and the maximum electron energy. Only converted cases 
were taken for averaging. 



Eq 


direction 


fraction of 


impart) 




{maxE^) 


(maxEe) 


[EeV] 




converted 






[EeV] 


[EeV] 


50 


strong B± 


137/1000 


263±189 


2.0±0.0 


6.8±3.2 


4.9±6.5 


100 


strong 


922/1000 


482±231 


2.01±0.11 


18±8 


2.2±5.3 




weak 


0/1000 


1 





100 





1000 


strong 


1000/1000 


3330±649 


9.7±2.3 


65±16 


4.7±7.1 




weak B±^ 


1000/1000 


544±143 


2.7±1.0 


167±66 


7.2±8.3 
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Table 3 



^max and RMS values for gamma induced showers of two different primary energies 
and arrival directions. Here the strong B± direction is defined as ^ = 53°, (f) = 177° 
and weak B_l as 6* = 53°, ^ = 357°. 



Eo [eV] 


direction 


fraction of con- 
verted 


{Xmax) [g/cm^] 


{RMS{Xmax)) 

[g/cm^] 


1019.5 


strong B± 


1/50 


1065 


90 


1020.0 


weak B± 


1/100 


1225 


175 




strong B± 


91/100 
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85 
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weak B± 


100/100 


1040 


40 




strong B± 


100/100 
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35 



E = 10^° eV 



E = 10^^ eV 
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CD 



-7 
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B=0.03 G 
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lg(hv/E) 
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Fig. 1. Bremsstrahlung spectral distribution: E = lO'^^ eV (left) and E = 10^^ eV 
(right) for different transverse magnetic field strengths: B = 0.3 G, 0.1 G and 0.03 
G for the uppermost, middle and lowest curves respectively. 



strong B: 9 = 60 deg, cp = 1 77 deg 
medium B: 9 = deg 
weak B: = 50 deg, (p = 357 deg 




100 1000 10000 



altitude a.s.l [km] 

Fig. 2. B± along different shower trajectories for the southern Auger Observatory. 
Curves from top to bottom: strong, medium, and weali field direction (see text). 



36 



X 
T3 



■D 



0.00035 
0.0003 
0.00025 



"B 0.0002 

O 
i_ 
Q. 







5e+19eV 






1e+20eV 


- / \ 




3.16e+20eV 






1e+21 eV 


/ 91% 


' 100% / 


, 100% \ 


12%y^ 







0.00015 

g 

I 0.0001 
\ 

5e-05 


2000 4000 6000 8000 10000 12000 
altitude a.s.l [km] 

Fig. 3. Altitudes of first conversion along strong field direction at Malargiie plotted 
for four different primary photon energies. The numbers attached to the curves 
indicate the total conversion probability (see text). 
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Fig. 4. Integrated conversion probability vs initial 7 energy for weak and strong field 
directions. 



37 



S 0.8 

!q 
ns 
J3 

2 0.6 

Q. 

I 0.4 
> 

8 0-2 



1 

S 0.8 
0.6 
0.4 
0.2 




50 EeV 



a 
n 
o 

Q. 

c 
o 
"co 



50 100 150 200 250 300 350 
azimuth [deg] 



100 EeV 



50 100 150 200 250 300 350 
azimuth [deg] 



1 

0.8 
0.6 
0.4 
0.2 


1 

£ 0.8 

CO 

2 0.6 

Q. 

I 0.4 

> 

§ 0.2 




1 . 1 1 , . 1 1 . , 1 1 1 1 1 1 1 






70 EeV 


r.-.r.sjJi'T, . . 1 . . . . 1 . . 


, 1 . , , . 1 , . .'ri^a.Tr.-c 



50 100 150 200 250 300 350 
azimuth [deg] 



1/ . - 




/ :','J 




: .";/ 

/ ;// 

'' ,'■■'/' 

■ 


150 EeV _ 


' . . . 1 . . . . 1 . . . 


\ K 

. 1 .... 1 .... 1 .... 1 .... r 



50 100 150 200 250 300 350 
azimuth [deg] 



Fig. 5. Total probability of 7 conversion for different arrival directions and four 
initial energies. Each curve corresponds to a different zenith angle: B = 80° for the 
uppermost curve down to 6 = 0° for the lowest one in steps of 10°. Computations are 
made for Malargiie magnetic conditions and azimuth 0° means the shower arriving 
from the geographic North. 
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Fig. 6. Energy distribution of photons (top) and electrons (bottom) in 1000 preshow- 
ers initiated by 10^° eV photons arriving at Malargiie from the strong field direction, 
compare Table 2. Spectra weighted by energy are plotted to the right. 
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Fig. 7. Distributions of Npart (top left), N^+^- (top right), max (bottom left) 
and max E^. (bottom right) for Eq = 10^^ eV and the strong arrival direction. 
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Fig. 8. Number of particles in the preshower for different altitudes of first 7 con- 
version. Plotted are the preshowers initiated by 10^° eV photons arriving from the 
strong B± direction. 
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Fig. 9. Energy carried by the preshower electrons at the top of atmosphere vs the 
altitude of the first 7 conversion for a primary photon energy of lO^'^ eV in the 
strong field direction. The three points in excess of the general trend are the cases 
where the primary photon converted at high altitude and one of the bremsstrahlung 
photons converted close to the top of the atmosphere. For further comments see the 
text. 
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